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An effective solution to the problem of Hermite interpolation with a clothoid curve is provided. At the beginning 
the problem is naturally formulated as a system of nonlinear equations with multiple solutions that is generally 
difficult to solve numerically. All the solutions of this nonlinear system are reduced to the computation of the 
zeros of a single nonlinear equation. A simple strategy, together with the use of a good and simple guess function, 
permits to solve the single nonlinear equation with a few iterations of the Newton-Raphson method. 

The computation of the clothoid curve requires the computation of Fresnel and Fresnel related integrals. Such 
integrals need asymptotic expansions near critical values to avoid loss of precision. This is necessary when, for 
example, the solution of interpolation problem is close to a straight line or an arc of circle. Moreover, some special 
recurrences are deduced for the efficient computation of asymptotic expansion. 

The reduction of the problem to a single nonlinear function in one variable and the use of asymptotic expansions 
make the solution algorithm fast and robust. 
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1. INTRODUCTION 



There are several curves proposed for Computer Aided Design | Farin(2002)[ |Baran 
et al.(2010)Baran, Lehtinen, and Popovic} De Boor(1978) |, for planning trajectories of 
robots and vehicles or the geometric design" of roads [De Cecco et al.(2007)De Cecco ' 



Bertolazzi, Miori, Oboe, and Baglivo; Scheuer and Fraichard(1997)^. The most important 



among these are clothoids also known as Euler's or Cornu's spirals, clothoi ds splines, i.e. 
a planar curve consisting in cl othoid seg ments, circle s and straight lin es, rDavis(1999)[ 
Meek and Walton (1992) : Meek and W alton(2QQ4) ; |Meek and Wal ton(2QQ9) ; Walton] 
and Mee k(2009)P , generalized clothoids or Bezier spirals pW alton and Meek(1996)]. 
Pjn:hagorean Hodograph [ Walton an d Meek(2007) ; [Farouki and Neff(1995)] , bi-arcs and 
conic curves are also widely used [[Pavlidis(l 983)17 It is well known that the best or most 
pleasing curve is the clothoid, despite its transcendental form. 

The procedure that allows a plane curve to interpolate two given points with as- 
signed (unit) tangent vectors is called Hermite interpolation, if given curvatures at 
the two poin ts is also satisfied, then this is called Hermite interpolation I jMcCrae and, 



Singh(2009) | . A single clothoid segment is not enough to ensure G^ Hermite interpolation, 
because of the insufficient degrees of freedom. Many authors have provided algorithms that 
use a pair of clothoids segments to reach the G^ Hermite interpolation requests. However, 
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often it is considered enough the cost-effectiveness of a Hermite interpolation, because 
it has been seen numerically that the discontinuity of the curvature is negligible. 
The purpose of this paper is to describe a new method for Hermite interpolation with 
a single clothoid segment, which do es not need to split the problem in mutually exclusive 
cases as in ^Meek and Walton(2009)| and does not suffer in case of degenerate Hermite data 
like straight lines or circles (see Figure[T]on the right). These are of course limiting cases but 
are treated naturally in our approach. Because of numerical stability in the computation, 
we preferred to introduce an appropriate threshold to avoid such degenerate situation in 
practice. A precise estimate of such switching points is discussed. Finally, the problem of 
the G^ Hermite interpolation is reduced to a solution of a single nonlinear equation, for ex- 
ample by the Newton-Raphson method. In order to provide a fast and accurate algorithm, 
we give a good initial guess for the Newton-Raphson method solver so that a few iterations 
suffice. 

The remainder of the article is structured as follows. In the next section we define the in- 
terpolation problem, in the third we describe the passages to reformulate it such that from 
three equations in three unknowns it reduces to one nonlinear equation in one unknown. 
This is enough from the theoretical point of view. The fourth section is devoted to the dis- 
cussion of a good starting point for the Newton-Raphson method, so that using that guess 
we achieve a quick convergence, and how to select a correct solution of the nonlinear equa- 
tion. Section [5] introduces the Fresnel integrals in their standard form and shows how the 
G^ interpolation problem (i.e. the argument of the trigonometric functions involved is not 
purely a second order monomial, but a complete second order polynomial) can be con- 
ducted to that form. Section [6] analizes what happens when the parameters of computation 
give numerical instabilities in the formulas, and ad hoc expressions for such cases are pro- 
vided. In the appendix there are the algorithms written in pseudo-code and a summary of 
the algorithm for the accurate computation of the clothoid spline. 



2. THE FITTING PROBLEM 

The construction of highway and railway routes and the trajectories of mobile robots can 
be split in the in construction of a piecewise clothoid curve which definition is given next. 

Definition 2.1 {Clothoid curve). The general parametric form of a clothoid spiral curve is 
the following 



x{s) ^ xq + J COS ^^^''''^ + KT + 1)0^ dr, 
vis) ^ yo + J sin Qw^'t^ + KT + -do^ dr, 



(1) 



where s is the arc length, n't -I- k is the linearly varjdng curvature, {xq, yo) is the starting 
point and i?o initial angle. Notice that ^k's^ -I- ks + z9o is the angle of the curve at arc length 
s. 



Clothoids curves a re computed via the Fresnel sine and cosine integrals (see Abramowitz 
and Stegun(1964) for a definition) which is discussed forward in Sectionjs] 



The determination of the parameters -do, k and k' are determined by points and angles at 
the extrema of the curve. Thus, the problem considered in this paper is stated next. 
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(3:0, J/0 ) 




Fig. 1. On the left Hermite interpolation schema. On the right some cases representing different final angles. 
The magenta arc is the limiting case of a circle. The other colored arcs represent intermediate cases. 

Problem 1. Given two points {xo,yo) cind {xi,yi) and two angles i?o and 'di, find a 
clothoid segment of the form ^ which satisfies: 



xo, 



x(0) 
x{L) = xi 



2/(0) = yo, 

y{L) = 2/1, 



arctan 
arctan 



x'(0) 
x'{L) 



^0, 



(2) 



with minimal positive L (the length of the curve). The general scheme is showed in Figure^on 
the left. 

Solution of Problem [T] is a zero of the following nonlinear system involving the imknowns 

L, K, k': 



F{L,k,k') 



/ xi — xq — COS (jk's^ + Ks + i9o) ds ' 
2/1 - 2/0 - Jq sin (^k's^ + ks + i?o) ds 
i?i - (i^'L^ + Ki + i9o) 



(3) 



Find the points such that F{L, k, k') = is a difficult problem to solve in this form. In the 
next section a reformulation of the problem in a simpler form permits to solve it easily. 

3. REFORMULATION OF THE PROBLEM 

Nonlinear system (|3]) is put in an equivalent form by introducing the parametrization s = tL 
so that the integrals involved in computation have fixed extrema independent of L: 

f B 2A\ ( cos {At^ + Bt + {}^) dr 

^[^'L' U) " ^V^Lfl sin {At^ + Bt + ^o) dr 

V i?i - (A + B + i?o) 

v^rhere A — ^k'L'^, B — Lk, Ax = xi — xq, Ay ~ yi — yo- 
The third equation in (|4j) is linear so that v^fe can solve it with respect to B, 



(4) 



B 



Ai9 ~A, Ai? = i?i - t9o, (5) 

and the solution of nonlinear system Q is reduced to the solution of the nonlinear system 
of two equations in two unknown, namely L and A: 



G{L,A) = 



Ax-L /q cos (At^ + (Ai) - A)t + i^o) dr ' 
Ay-L /g sin {At^ + (Ai? - A)t + ^0) dr ^ 



(6) 
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followed by the computation of B using (js). Nonlinear system ([6j) is easier to solve than Q. 
We can perform further simplification using polar coordinates to represent (Ax, Ay), 
namely 

Ax = r cos 93, Ay = r simp, (7) 

and from ([7]) and L > we can define two new nonlinear functions f{L, A) and g{A), 
where g{A) is independent of L, as follows: 

Using the identity sin(Q! ~ /3) = sin a cos /3 — cos a sin /?, function g{A) simplifies in: 

g{A) = e{A;Ad,A^), (8) 

where Atp — 'do — ip and 

e(A; Ai9, Av5) = /" sin (At^ + (A?9 - + A(^) dr, (9) 

while using the identity cos(a — /3) = cos a cos /3 + sin a sin /3, function /(L, ^) reduces to 

f{L,A) = ^Ax^ + Ay^~Lh{A), 

h{A) = cos (Ar^ + {Ad - A)t + A(/3) dr = 6 (A; Az?, Ayp + |) . 

Remark 3.1. Defining 0o = ^^o — and 01 = -iJi — g p we have A -d = 0i — (^0 and A^p — <po, 
moreover, the angles 0o and are the same used in Walton and Meek(2Q09), to derive 
interpolant. 

Lemma 3.2. The function 9 defined in ^ satisjy 

Q{A;Ad,z) =0, ^ e (A; Ai?, z + I) 7^ 0. (11) 



Proof. If A = by a simple computation we obtain 

cos A?9 - cos(Ai9 + z) 

e(A;Ai?,z) = ) 

^ ' ' 1 — cos z 



if At? 7^ 0, 
if Ad = 0, 



and it is immediate to verify dll] ) . Let A > 0, after some manipulation we have 

V2A, 

where 



-<d{A-Ad,z) = (5(a) +5(6)) cos ?7- (C(a) + C(6)) sin r/ 



^-Ai? , A + At9 TT 2 

a= — , 0= — , , 7] ~ —a — z. 

V2TrA 2 

If implication ([TT]) is false, i.e. both 9 {A; A-d, z) = and 9 {A; Ad, z + tt/2) = then from 
previous equations it follows 

S{a) + S(b) = 0, C{a) + C{h) ^ 0, 

and the symmetry C(—z) = — C{z) With. S{—z) = —S{z) implies that points {C{a),S{a)) and 
{C{-b),S{-b)) are coincident points on the Comu spiral in parametric form (C(t), S{t)) and 
thus a~—b. Hence, 

, A -Ad A + Ad V2A 

= a + 6= H ; = 

V27rA V27rA Vtt 
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and thus A = which contradict the assumption that A > 0. The case with A < is 
similar. □ 

Actually, the situation can be understood graphically looking at the plots oi g{A) and h{A), 
as in figure |2]on the left:. It is possible to reverse the implication by the following lemma. 

Lemma 3.3. The solutions of the nonUnear system (|6]) are given by 

Ay^ Ad -A , 2A 



L = 



L 



(12) 



Proof. Let L, A satisfy ^ and ([12]). Then f{L,A) = and thus G{L,A) = 0. From 
Lemma 3.2 when g{A) = then h{A) =^ and thus L is well defined. □ 



where A is a root of g{A) defined in equation ^ and h{A) is defined in ( [IQ] ). 



Remark 3.4. The solution of problem [T] is build using Lemma 3.3 selecting solution of 
minimal length. Notice that the parameter A corresponding to minimal (positive) L depend 
only on the relative angles and it is independent on scaling factor So that we can compute 
A using only A^p and A-d assuming r 



Ay2 ^ 1. 



Hence the interpolation problem is reduced to a single (nonlinear) equation that can be 
solved numerically with Newton-Raphson Method. By a graphical inspection we notice 
that the roo ts o f g{A) are si mple, so that Newton-Raphson converges quadratically. Func- 

implements a simple strategy for the computation of minimum 



tion f indA 



with 



buildGrid 



length clothoid on agrid oif points and angles. This algorithm is used to compute surface 
presented in Figure [3]on the top. 



Function findA(^guoss, A-d, Aip, tol) 



1 Define g as g{A) := 6 (A, Ai?, A(/?); Set A ^gucss; 

2 while > tol do A A - g{A) / g' {A) ; 

3 return A; 



Function buildGrid (iV, Af ) 


1 forj 


= 0, 1, . . . , A/ do Af^^^ ^ 07Af)40 - 20; 


2 for i,j — 


0, . . . , iV do 


3 




-2ni/N-TT; A'd ^ 2nj /N - n; Lij oo; 


4 


for k 


= 1, . . . , A/ do 


5 




if QiAf"^', Ai9, Avp) e(yl|rr, Ai9, A^) < then 


5 






A ^ f indA ((Af^'^" + AlTi")/2,A^A'p); 


7 






L ^ l/e(I, Ai9, Avp + 7r/2); 


8 






if L > and L < Lij then Lij L; Aij A; 


9 




end if 


10 


end for 


11 end for 




12 return A, L; 



In the next section a recipe to select a good initial point is described. 

4. SOLUTION OF HERMITE INTERPOLATION PROBLEM 

By a graphical inspection the so lution is in the range [—20, 20] for angles Aip and At? in 
the interval [— vr, tt]. The function buildGrid|is used to compute all the possible solutions in 



this interval for angles ranging in [—tt, ttJ and selecting the ones with minimum length. The 
strategy is very simple; the interval [—20, 20] is divided in small intervals and where g{A) 
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y 

li 






mm 




Fig. 2. On the left: plot of g (red line) and h (bl ue l ine), when g{A) = then h{A) 7^ so that the length L is 
well defined. In centre clothoids computed using \13\ as starting point for Newton-Raphson iterative scheme, on 
the right the minimum length clothoids. The clothoids start from (0, 0) angle to (—0.95, 0.31) and angle ranging 
from — TT to IT. 



changes sign Newton-Raphson iterative scheme is used. If convergence is attained and L, 
computed using ( |12D , is positive, this solution is a candidate to be the final solution. 

This strategy, although simple, is effective and permits to evalu ate the solut ion of Prob- 
lem[T]at any point. Figure [sjshov^fs the results of computation using |buildGrid v^ith iV = 32 
and M ~ 128. Notice that the solution is discontinuous for higher and lower angles. To 
avoid this jump in the solution we drop the requirement of minimal length for all angles 
and solve the following problem. 



Problem 2. Given two points {xo,yo) cind(xi,yi) and two angles do cind -di, find a 
clothoid segment of the form ([T]) which satisfies ^ with minimal positive Lfor |i?o| and 
less then tt/2. For greater angles we choose the solution that make the function A{Ai}, Aip) 
with A{A-d, Aip) e {A \ Q{A, Ad, A(p) = 0} a continuos surface. 

Looking at Figure |3] we guess that Ai^Ad, A(p) can be approximated by a plane. A simple 
fitting on the solution of Problem |2] with a plane, results in the following approximation for 

A{Ad,Aip): 



A{Ai), Aip) « 2.4674 Ad + 5.2478 Aip. 



(13) 



Using ( |13D as starting point for Newton-Raphson the correct solution for Problem |2]is found 
in few iterations. Figure[3]show the computed solution found using this strategy. Computing 
the solution with Newton-Raphson starting with the guess given by ( \13} in a 1024 x 1024 
grid with a tolerance of 10~^° results in the following iteration distribution: 



iter 


1 2 


3 


4 


5 


#cases 


1 212 


148880 


846868 


54664 




less than 1% 


14% 


80% 


5% 



Notice that the average number of iteration is 3.9 while the maximum number of iterations 
per point is 5. Figure [2] on the right shows the difference of computing minimum length 
clothoid curve and clothoid curve of Problem |2] The minimum length produces a bad look- 
ing solution for some angles combinations. The curve produced by the solution of Problem 2] 
looks better and in any case curve produced by the solution of Problem [T] and Problem 2] 
are identical for a large combination of points and angles. 
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T he algorithm for th e clothoid computation is resumed in function buildClothoid which 
uses normalizeAngle to put angles in the correct ranges. 



Function normalizeAngle(</j) 


1 while if > +7r do (/3 <^ 


- ip-2n; 


2 while ip < —IT do (/? -f 


-<p + 2n; 


3 return ip; 





the equation g{A) 



by calling 



Function |buildClothoid| solves eq uation ([8]) that is it finds a solution of 

described at the end of Section [s] 



f indA 



Function buildClothoid (xq, yo, "^o, xi, yi, 

1 At xi - xo; Ay <— yi — yo; 

2 Compute r and (p from r cos ip = Ax and r sin p — Ay; 

3 Ap-(r- normalizeAngle('!9o — p); A'd-i— normalizeAngle(i?i — i9o); 

4 A ^ findA(2.4674 Ai? + 5.2478 A(/p,Ai9,A(/p); 

5 L^r/e{A,A-d,Ap+^); n ^ [A-d - A) / L; k'^{2A)/L^; 

6 return /t, k, L 



The algorithm needs accurate computation of Fresnel related functions g{A) and h{A) 
with relative derivative which are discussed in the next section. 





Fig. 3. On the top, computation of parameter A and minimal length and for nonlinear system j6]l. On the bottom 
comput ation of para meters A and L for Problem[2] Computation of the figures on the top are done using auxiliary 
function lbuildGridl 
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5. COMPUTATION OF FRESNEL AND RELATED INTEGRALS 

Among th e various possible definitions, w e choose the following one, which is the same 



adopted in Abramowitz and Stegun(1964) 



Definition 5.1 (Fresnel Integrals). Following Abramowitz and Stegun(1964) the Fresnel 
Integrals are defined as 

C{t) = cos (|t2) dr, S{t) = sin (^r^) dr. (14) 
Remark 5.2. Some authors prefer the definition 



Cit) = 



f cos{T^)dT, and 5(t) = / sm{T^)dT. 
Jo Jo 



However by using the identities 



it is easy to pass between the two definitions. 



To compute the stand ard Fresnel integrals o ne can use algorithms described in Sny- 
|der(1993) , Smith(2011) and Thompson(1997)' or using continued fraction expansion as 
in Backeljauw and Cuyt(2009) It is possible to reduce the integrals ([T]) to a linear combina- 
tion of this standard Fresnel integrals ( |141 ) . However simpler expression are obtained using 
also the momenta of the Fresnel integrals: 

Ck{t) = jy cos (|t2) dr, Skit) - jy sin (^r^) dr. (15) 

Notice that C{t) :— Co{t) and S{t) :— So{t). Closed forms via the exponential integral or 
the Gamma function are also possible, however we prefer to express them as a recurrence. 
Integrating by parts, the following recurrence is obtained: 



Ck+i{t) - - (t'^sin^^t^'^ -kSk-i{t) 
Sk+i{t) - ^ (fcCfc_i(t) - t'' cos 



(16) 



Recurrence is started by computing standard Fresnel integrals ( |14D and (by using the 
change of variable z = r^) the following values 

Ciit) - isin (|t2) , Slit) = l{l- cos , 

From recurrence follows that Ck{t) and Skit) with k odd do not contain Fresnel integrals 
( |14| ) and are combination of elementary functions. 

The computation of clothoids relies most on the evaluation of integrals of kind ([9]) with 
their derivatives. The reduction is possible via a change of variable and integration by parts. 
It is enough to consider two integrals, that cover all possible cases: 

Xkia,b,c) ^ / r''cos(^r^ + 6r + c) dr, 

^ a ^'^^ 

Yfc(a, fe,c) J r*" sin (^^r^ + 6r + c^ dr. 

In fact, integrals ([s]) and ( |1QD with their derivatives can be evaluated by using the identity 
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where a = k' s^, b — ks, and c = ??o so that 

gir,) = Yo{2r,, - 77, ~ip), hir,) = Xo{2r,, - r,, -^), 
and finally, equation can be evaluated as 

x{s) = Xq + s Xo{k' s'^ , KS,do), 

y{s) = yo + sYo{K's'^,Ks,^o)■ 
FTom now on, algorithms for accurate computation of X^. and are discussed. From the 
well known identities 

cos(a + (3) — cos a cos f3 — sin a sin /?, sin(a + /?) = sin a cos f3 + cos a sin /3, (18) 
integrals ( [17D can be rewritten as 

Xfe (a, &, c) = Xk (a, 6, 0) cos c — Yfe (a, 6, 0) sin c, 
Yfe (a, b, c) = (a, 6, 0) sin c + Yk{a,b, 0) cos c, 

Thus defining Xk{a,b) := Xfe(a, 6, 0) and Yk{a,b) :— Yfc(a, 6,0) the computation of ( |171 ) 
is reduced to the computation of Xk{a,b) and Yk{a,b). It is convenient to introduce the 
following quantities 

so that it is possible to rewrite the argument of the trigonometric functions of Xk{a, b) and 
Yfc(a,6) as 

a o , TT / ■v/|a| sb \ sb'^ ix , ^.2 
2 2 \^ V7r|a|y 2 |a| 2 ' 

By using the change of variable ^ ^ t z + i with inverse r = z^^ — £) for (a, 6) and the 
identity ( |18P we have: 

Xk{a,b) = J z-''X^-^)^cos(^e'+7) dC, 



z 

A; 



3=0 



where 



ACj ^ C,{£ + z) - Cj{£), ASj = Sj{e + z) ~ Sj(£), (19) 

are the evaluation of the momenta of the Fresnel integrals as defined in ( |15D . Analogously 
for Yk{a, b) we have: 

Ykia, b) = z-''-^ J2 (^^^''^^ [sin7ACj + scos7A5j] . 

A cheaper way to compute Xk{a, b) and Yk{a, b) is via recurrence as in ( |16| ) for the compu- 
tation of Fresnel momenta. Integrating the identity, 

A sin (^^2 + = fci^-i sin (^t^ + 6t) + + 6) cos (^t^ ^ 
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we have 



3in + = fcYfc_i(a, b) + aXk+i{a, b) + bXk{a, b) 



which can be solved for Xk+i{a,b). Repeating the same argument with cos((a/2)i^ + bt) 
and solving for ¥^+1(0, h) we obtain recurrences for X and Y: 



Xk+i{a-b) 



- sin - + 



Yk+i{a.b) = i (fcXfc_i(a,6)-&rfc(a,&)-cos(^+6)) . 



(20) 



From the identities 



/ (ar + b) cos {^t"^ + brj dr aXi{a, b) + 6Xo(a, b) = sin {^+bj , 
J {ar + b) sin (^%^ +bTj dr ^ aYi{a,b) + bYo{a,b) =l-cos(^^+6^ 



(21) 



(22) 



it follows 

Xi(a,6) = ^ (sin(^+&) -feXo(a,6)) , 
Y,{a, 6) = i (1 - cos + fe) - bYoia, 6)) . 
To complete recurrence we need the initial values Xo{a, b), Ya{a, b): 

Xo{a,b) — z~"^(cos7ACo — ssin7AiSo), 
Yo{a,b) = z^"^ ( sin 7AC0 + sees 7A5o), 

where Co and So is defined in ( [191 ). This recurrence may be inaccurate when |a| is small, in 
fact z appears in the denominator of many fractions. For this reason for small values of \a\ 
we substitute this recurrence with asymptotic expansions. 

6. ACCURATE COMPUTATION WITH SMALL PARAMETERS 

When the parameter a is small we use identity ( |18| ) to derive series expansion 

^1 

Xkia,b) 



r'" cos ( -T^ + 6t ) dr 



fa 



dr 



E 



(-1)" /a\2n 



- {2n)\ 



Tl = 

00 



1-^ I— ir° /a\2n+l 
^4„+.(0,6) - 5: ^^-^ (-) r4„+2+fc(0,6) 
ri=0 ^ ^ ''• 
aKln+2+fc(0, fo) 



V (-1)" 

(2n)! V2y 



(23) 



X 



in+k 



(0,&) 



2(2n + 1) 



and analogously using again identity ( |18D we have the series expansion 

»i 

Yfc(a,6) = / r*^ sin f "r^ + 6r ) dr 



"(I 

(-1)" /a\2" 



(24) 



^ (2n)! 



Y, 



4n+k 



(0,6) 



aX4„+2+fc(0, 6) 
2(2n + l) 
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From the inequalities 



1 



fc + 1 



k + 1' 



we can estimate the remainder for the series of Xk 

(-1)" /a\2« 



E 



(2n)! V2 
- ^ (2nV. V2/ 



n—p 
oo 



n—p 



^4„+fe(0,6) 



1 



2(2n+ 1) 



- (2) E (2{n-p))\ (2 



4n + l 2(2n + l)(4n + 3) 
a\ 2(n-p) 



^(2) Et^G) -(2) 



n—p 

00 



1 / a \ 2" / a \ 2p 



(2n)! V2 

n=0 ^ ^ 

The same estimate is obtained for the series of Yk 



(25) 



Remark 6.1. Series ( |23| ) and ( |241 ) converge fast. For example if \a\ < 10 and p — 2 
the error is less than 6.26 • 10^^^ while if p = 3 the error is less than 1.6 • 10^^^. 

Recurrence (|2Q])-(|2T])-(|22]) permits to compute Xk{a,b) and Yk{a,b) at arbitrary precision 
when a 7^ 0, but when a = it modifies to 

Xfc(0,fe) = b-\sinb-kYk-i{G,b)), 
Yu{0,b) = 6-i(fcXfc_i(0,6)-cos6), 
with starting point 

Xo(0,6) = 6-1 sin 6, Yo{Q.b) = b'^ {I - cosb) . 

Recurrence ( |25] ) with ( |23D and ( [241 ) permits computation when 7^ 0, unfortunately recur- 
rence ( |25] ) is highly unstable and cannot be used. In alternative an explicit formula based 
on Lommel function s^^y{z) can be used | ]Shirley and Chang(2003)p . Explicit formula is the 
following 



Xfc(0,6) 



cos 



(1 + k)b''+^ 



1 + fc' 



, , fcSfc+3 J (6) sin& + g(&)sfc+i^ 1 (b) sin 6 



(26) 



(2 + fc)6'=+3 

where = b~^ sin 6 — cos 6 and g{b) = f{b){2 + fc). Lommel function has the following 
expansion (see http : //dlmf . nist . gov/11 .9^ 

(-z2)" 



^(m,^)= n(('^ + 2m-l)2-z.2)^ (27) 



and using this expansion in ( |26| ) results in the following explicit formula 

Xfe(0,6) = . 3(6) + B(6)u,,+ 3 .(6)+ 



yfc(0,5) = C(6K+3 3(6)+i?(6K+i 1(6) + 



1 + fc' 
sin 5 

2 + fc' 
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where 



2\n 



C{b) 



Aib) = 



kb sin b 

1 + k ' 



Bib) 



(sin 6 — & cos 6)6 



1 + fc 



6^ sin 6 
2 + fc ' 



D(b) = sin 6—6 cos 6. 



7. CONCLUDING REMARKS 

The proposed algorithm is detailed in the appendix using pseudocode and can be easily 
translated in any programming language. A MATLAB implementation is furnished to test the 
functions introduced in the paper It is important that the computation of Fresnel integrals 
is accurate because they are involved in many expansions for the clothoid fitting. In the 
MATLAB implementation of the proposed algorithm for Fresnel integrals approximation a 
slightly modified version of Venkata Sivakanth Telasula script, available in MATLAB Central, 
was used. 



A. ALGORITHMS FOR THE COMPUTATION OF FRESNEL RELATED INTEGRALS 

We present here the algorithmic version of the analytical expression we derived in Sec- 
tion [5] and [6} T his algorithm are necessary for the computation of the main function 



buildClothoid 



^of Section [4] which takes the input data (a;o, yo, -do, x\, yi, di) and returns 

the para meters (k, k' , L) that solve the problem as expressed in equation ([l]). Function 
evalXY computes the generalized Fresnel integrals ( |171 ) . It distinguishes the cases of a larger 
or s maller than a threshold e. The value of e is discussed in Section [6j see for example Re 
Formulas ([20]) -([21]) -([22]), used to com pute Xk (a 
> e, are implemented in function 



mark 6.1 
cision when 



evalXY aLarge 



to comp ute Xkia,b) and Yk{a, 6) at arbitrary precision when 



function 



evalXYaSmall 



in function 



b) an d Yk{a, 6) at arbitrary pre- 
Formulas ([23) -([24]), used 
< £, are implemented in 



evalXYaZero 



function IrLommelF 



This function requires computation of ( |25| ) a nd (|26p implemented 
which needs (reduced) Lommel function ([27^ implemented in 



Function evalXY(a, 6, c, fc) 



1 if|a| <ethen ^- [evalXYaSmall| :a,b,fc,p); else X°,y" ^- [evalXYaLarge| :a,b,fc); 

2 for j = 0, 1, . . . , fc do 

3 I Xj cos c — Yj' sin c; Yj ^ X° sin c 4- cos c; 

4 end for 

5 return X, Y 
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Function evalXYaLarge(a, b, fc) 



a -v/tt „ sb sb 1 , 

ei ; Ji ;r-p-r; t ^ -a + b; 



\a\ ^''^ '^\'^\ 2 

2 ACo C(t + z) - Cil); 
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Function rLommel((U, i^, 6) 



It (/i + + 1) """(/i — + 1) r ^ t; ri <— 1; 

(-fa) fa 

'2n + ^- i/ + 12n + + + 



J (-fa) fa 
2 while |t| > e do t ^ t- ^ — -; r^r + t; + 



3 return r 



Function evalXYaZero(6, fc) 



1 if |b| < e then 

b^ f 6V faV fa^ 

2 " " 

3 else 

4 



sin 6 , , 1 — cos b 
Xoi 7—; Yo ■ 



b ' " b • 

5 end if 

6 ^ fosinfe; D -(r- sinb — 6cos&; B <— &_D; C — fe^sinfc; 

7 for fc = 0, 1, . . . , fc — 1 do 



fcA lrLommellf fc + |. §. + g ELommell ffc + ^ 

Xk+ii ^ ^ ^ 



2 , 2 , fa I + cos b 



1 + k 

3 3 



C IrLommell f fc + § . I . b J + sin b / ^ ^ 

Yk+ii 2 + k +DlrLominellf fc+ - . -,6); 

10 end for 

11 return X, Y 
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Function eva]XYaSmall(a, b, k, p) 



1 X°,Y° 

2 for j = 0,1 

3 for n = 1 , 2 



evalXYaZero 



b,k + 4p + 2); t^l; 



, k do Xj 
,pdo 



Y 



j+2> 



16n(2n- 1)' 
for ji = 0, 1, . . . , fc do 



end for 

8 end for 

9 return X, Y 



4n+j in+j + 2 . 

4n + 2 ' 



yO I • 



4n+j 



4n + 2 



An+j + 2 . 
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